{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python Practice 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Practicing ANOVA on Random Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's play a bit with the ANOVA analysis on random data. First, we will generate a random dataset of 50000 points and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_points = # put your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the random points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points indeed look pretty random. Since ANOVA test is designed for multiple data samples (or multiple groups), we will group the data into 50 samples of 1000 points each. \n",
    "Next, we will calculate the mean and standard deviation for each sample and plot them with error bars.\n",
    "\n",
    "Instead of creating the new array with the size of (50, 1000) we can use the `np.reshape` function to reshape the original array into the desired shape.\n",
    "It works in the following way:\n",
    "\n",
    "```python\n",
    "N = 1000\n",
    "M = 50\n",
    "initial_array = np.random.normal(0, 1, N*M)\n",
    "reshaped_array = np.reshape(initial_array, (M, N))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 50\n",
    "num_points_per_sample = 1000\n",
    "random_samples = # put your code here\n",
    "\n",
    "print(\"Old shape: \", random_points.shape)\n",
    "print(\"New shape: \", random_samples.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = # put your code here\n",
    "stds = # put your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the means values of each sample with the errorbars defined by their standard deviation, we can use the `plt.errorbar` function.\n",
    "\n",
    "```python\n",
    "plt.errorbar(x, y, yerr=std, fmt='o')\n",
    "```\n",
    "\n",
    "Where `x` is the x-axis values, `y` is the y-axis values, `yerr` is the error values, and `fmt` defines the symbol for plotting the points (in this case, 'o' for circles).\n",
    "\n",
    "Since in our case we don't have specific x-axis values, we can use the indices of the samples as x-axis values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the means and standard deviations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we expected for the standard normal distribution, the mean values are around 0, and the error bars are about 1.0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running ANOVA test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now perform the ANOVA test on the random data. The ANOVA test is used to determine whether there are any statistically significant differences between the means of three or more independent (unrelated) groups.\n",
    "\n",
    "Given that we have 50 samples ($n$) of 1000 points each ($m$), we need to calculate the following values to perform the ANOVA test:\n",
    "\n",
    "\n",
    "- The error sum of squares (SSE): $$SSE = \\sum_{i=1}^{n} \\sum_{j=1}^{m} (X_{ij} - \\bar{X}_{i})^2 $$ \n",
    "We can calculate the SSE by first calculating the difference between the data in each sample and the mean of that sample, squaring the difference, and then summing all the squared differences.\n",
    "\n",
    "- The sum of squares between the samples (SSB):  $$ SSB = \\sum_{i=1}^{n} m (\\bar{X}_{i} - \\bar{X})^2 = m \\sum_{i=1}^{n}  (\\bar{X}_{i} - \\bar{X})^2 $$\n",
    "We can calculate the SSB by first calculating the difference between the mean of each sample and the overall mean, squaring the difference, multiplying by the number of data points in each sample, and then summing all the squared differences.\n",
    "\n",
    "- The total sum of squares (SST): $$ SST = SSE + SSB $$\n",
    "\n",
    "- The degrees of freedom between samples: $$ df_{b} = n - 1 $$\n",
    "\n",
    "- The degrees of freedom within samples: $$ df_{e} = n(m-1) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement these functions:\n",
    "\n",
    "def calculate_sse(samples, means, num_samples):\n",
    "    sse = # put your code here\n",
    "    return sse\n",
    "\n",
    "def calculate_ssb(means, num_samples, num_points_per_sample):\n",
    "    ssb = # put your code here\n",
    "    return ssb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SSE = calculate_sse(random_samples, means)\n",
    "SSB = calculate_ssb(random_samples, means, num_points_per_sample)\n",
    "\n",
    "dfb = # put your code here\n",
    "dfe = # put your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Error Sum of Squares (SSE): \", SSE)\n",
    "print(\"Sum of Squares between samples (SSB): \", SSB)\n",
    "print(\"Degrees of freedom between samples (dfb): \", dfb)\n",
    "print(\"Degrees of freedom within samples (dfe): \", dfe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the vectorized operations in NumPy to calculate these values efficiently. The main idea is to directly operate on 2D array (`samples`) to calculate the mean values and then use broadcasting to calculate the differences and sums.\n",
    "\n",
    "The issue is that `means` is a 1D array, and we cannot directly subtract it from the 2D array `samples`. We first need to reshape the `means` array to have the same shape as the `samples` array.\n",
    "\n",
    "```python\n",
    "means = means.reshape(len(samples), 1) # Reshape to (n, 1)\n",
    "```\n",
    "\n",
    "Now we can subtract the `means` array from the `samples` array. It will automatically subtract the corresponding mean value from each row of the `samples` array.\n",
    "\n",
    "```python\n",
    "diff = samples - means\n",
    "```\n",
    "\n",
    "This will return the array of the same shape as the `samples` array, where each element is the difference between the data point and the mean of the corresponding sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_sse_vectorized(samples, means):\n",
    "    return np.sum((samples - means.reshape(len(samples), 1)) ** 2)\n",
    "\n",
    "def calculate_ssb_vectorized(means, num_points_per_sample):\n",
    "    total_mean = np.mean(means)\n",
    "    return num_points_per_sample * np.sum((means - total_mean) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SSE_vec = calculate_sse_vectorized(random_samples, means)\n",
    "SSB_vec = calculate_ssb_vectorized(means, num_points_per_sample)\n",
    "\n",
    "print(\"SSE and SSE_vec are equal: \", np.isclose(SSE, SSE_vec))\n",
    "print(\"SSB and SSB_vec are equal: \", np.isclose(SSB, SSB_vec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To finish the ANOVA test, we need to calculate the mean square between samples (MSB) and the mean square within samples (MSE). The F-statistic is then calculated as the ratio of MSB to MSW.\n",
    "\n",
    "$$ MSB = \\frac{SSB}{df_{b}} $$\n",
    "$$ MSE = \\frac{SSE}{df_{e}} $$\n",
    "$$ F = \\frac{MSB}{MSE} $$\n",
    "\n",
    "Finally, we compare the F-statistic to the critical value from the F-distribution with degrees of freedom $df_{b}$ and $df_{e}$ at a desired significance level (in our case 0.05)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MSB = # put your code here\n",
    "MSE = # put your code here\n",
    "F = # put your code here\n",
    "\n",
    "print(\"Mean Square between samples (MSB): \", MSB)\n",
    "print(\"Mean Square within samples (MSE): \", MSE)\n",
    "print(\"F-value: \", F)\n",
    "\n",
    "print(\"Critical F-value: \", 1.354) # for alpha = 0.05, dfb = 49, dfe = 49950"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our F-statistic is less than the crititical value, we cannot reject the null hypothesis (i.e. we cannot say that the data in the samples differs significantly).\n",
    "\n",
    "Let's summarize the results in the [ANOVA table](https://d138zd1ktt9iqe.cloudfront.net/media/seo_landing_files/anova-table-1642579664.png)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_anova_table(SSE, SSB, dfb, dfe):\n",
    "    MSB = SSB / dfb\n",
    "    MSE = SSE / dfe\n",
    "    F = MSB / MSE\n",
    "    total_SS = SSB + SSE\n",
    "    total_df = dfb + dfe\n",
    "\n",
    "    print(f\"{'Source':<20} | {'Sum of Squares (SS)':<20} | {'Degrees of Freedom (df)':<20} | {'Mean Square (MS)':<20} | {'F-value':<10}\")\n",
    "    print(\"-\" * 105)\n",
    "    print(f\"{'Between samples':<20} | {SSB:<20.2f} | {dfb:<20}    | {MSB:<20.2f} | {F:<10.2f}\")\n",
    "    print(f\"{'Within samples':<20} | {SSE:<20.2f} | {dfe:<20}    | {MSE:<20.2f}\")\n",
    "    print(f\"{'Total':<20} | {total_SS:<20.2f} | {total_df:<20}\")\n",
    "\n",
    "\n",
    "print_anova_table(SSE, SSB, dfb, dfe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Practicing ANOVA on Lab Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's repeat the same analysis on the data from our lab. This data will contain some slow drifts, which might introduce some differences between the samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}